function vmd_potential (fbase, fnum)
%
%   vmd_potential(fbase, fnum)
%
%   Reads density file. Writes potential and electric field.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Parameters for background potential and system size are edited in the file.


%---------------------------------------------------------

%-- system size

jx=48; jy=48; jz=16;

%-- physics paramters

Lx=450; 
Ly=450; 

phi00        = 0.01; 
SpeckleWidth = 99.91139; 
SpeckleDepth = SpeckleWidth;
kz           = 0.34; 
nz           = 1;


%-- colormap bounds

maxEz = 0.05;  maxEx = 0.005; maxEy = 0.005; maxPhi=0.1; maxPh0=0.01


%---------------------------------------------------------


%-- set background potential --

  dx = Lx/jx;
  dy = Ly/jy;
  dz = 2*pi*nz/(jz*kz);

  x  = (0:jx-1)*dx - Lx/2;
  y  = (0:jy-1)*dy - Ly/2;
  z  = (0:jz-1)*dz;

  [xx,yy]= meshgrid(x, y);

  rmod = exp( -20*(xx/SpeckleWidth).**2  - 20*(yy/SpeckleDepth).**2 );

  phi0=zeros(jz,jy,jx); 
  for iz=1:jz
    phi0(iz,:,:) = rmod * phi00 *exp(i*kz*z(iz)) ;
  end

  phi0FT = fftn(phi0);


%-- read density distribution

  fn = num2str(fnum,'%04d');

  fname = [fbase, '.dn.', fn];

  fid = fopen(fname, 'rb');
  rho = fread(fid, jx*jy*jz, 'double');
  fclose(fid);

  rho = reshape(rho, jz,jy,jx);


%-- solve poisson equation 

  rho = fftn(rho);

  dkz = kz/nz;
  dkx = 2*pi / Lx;
  dky = 2*pi / Ly;

  kx = [0:jx/2,-jx/2+1:-1] *dkx;
  ky = [0:jy/2,-jy/2+1:-1] *dky;
  kz = [0:jz/2,-jz/2+1:-1] *dkz;

  [kky, kkz, kkx] = meshgrid(ky, kz, kx);

  kk = kkx.*kkx + kky.*kky + kkz.*kkz ;

  phi = rho./kk;

  phi(1,:,:) = 0;
  phi = phi + phi0FT;
  phi(:,1,1) = 0;            % spatially uniform in x not allowed (03/21/2011)


  Ex  = -i*kkx.*phi;  
  Ey  = -i*kky.*phi;  
  Ez  = -i*kkz.*phi;  

  Ex  = real(ifftn(Ex));
  Ey  = real(ifftn(Ey));
  Ez  = real(ifftn(Ez));
  phi = real(ifftn(phi));
  phi0 = real(phi0);


%-- images  --------------------------------

  fname = [fbase, '.Ex.', fn];
  img_slices(Ex, fname, -maxEx, maxEx)

  fname = [fbase, '.Ey.', fn];
  img_slices(Ey, fname, -maxEy, maxEy)

  fname = [fbase, '.Ez.', fn];
  img_slices(Ez, fname, -maxEz, maxEz)

  fname = [fbase, '.phi.', fn];
  img_slices(phi, fname, -maxPhi, maxPhi)

  fname = [fbase, '.ph0.', fn];
  img_slices(phi0, fname, -maxPh0, maxPh0)

end


%---------------------------------------------------------
